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We investigate thermodynamic properties like specific heat cy and susceptibility y in anisotropic 
J 1 -J 2 triangular quantum spin systems {S = 1/2). As a universal tool we apply the finite tem¬ 
perature Lanczos method (FTLM) based on exact diagonalization of finite clusters with periodic 
boundary conditions. We use clusters up to A = 28 sites where the thermodynamic limit behavior 
is already stably reproduced. As a reference we also present the full diagonalization of a small 
eight-site cluster. After introducing model and method we discuss our main results on CyiT) and 
X(T)- We show the variation of peak position and peak height of these quantities as function of 
control parameter J 2 /J 1 . We demonstrate that maximum peak positions and heights in Neel phase 
and spiral phases are strongly asymmetric, much more than in the square lattice J 1 -J 2 model. Our 
results also suggest a tendency to a second side maximum or shoulder formation at lower tempera¬ 
ture for certain ranges of the control parameter. We hnally explicitly determine the exchange model 
of the prominent triangular magnets CS 2 CUCI 4 and Cs 2 CuBr 4 from our FTLM results. 


I. INTRODUCTION 

The 2D triangular S' = 1/2 Heisenberg antiferromag- 
net (HAF) is the archetypical geometrically frustrated 
quantum magnet [T]. Therefore it has been studied 
in countless investigations. Although there is no pure 
physical realization due to in-plane symmetry breaking 
and inter-plane coupling several compounds fall into the 
wider class of anisotropic triangular magnets where the 
exchange coupling J 2 along one of the triangle sides is dif¬ 
ferent from Ji along the other two sides [2]. In fact this 
generalization is not a drawback but opens interesting 
possibilities. It allows to consider a smooth interpolation 
between square lattice HAF to isotropic triangular mag¬ 
net to quasi-lD chain compounds as function of a single 
tuning parameter J^jJ\- 

Most of the work on this magnetic system was focused 
on the zero temperature phase diagram as function of 
anisotropy ratio, spin wave excitations and influence of 
quantum fluctuations as well as effects of external helds 
like plateau formation in the magnetization. We will not 
touch these fundamental topics but refer to Ref. [5] and 
references cited therein. Here we will focus on a more 
practical issue, namely the determination of the exchange 
constant or at least their anisotropy ratio from thermody¬ 
namic quantities of triangular magnets. This can be done 
by analyzing the experimental temperature dependence 
of specific heat and magnetic susceptibility where the 
former is less favorable due to lattice contributions [3]. 
A simplified but not unambiguous determination of ex¬ 
change constants is possible by identifying position and 
height of the maximum in thermodynamic quantities and 
comparing with theoretical predictions. Further useful 
tools of diagnosing the exchange model are magnetother- 
mal, e.g. magnetocaloric properties as well as high field 
magnetization and saturation field. This program has 
sofar mainly been carried out for frustrated J 1 -J 2 square 
lattice magnets [IHS]. 

Here we focus on the numerical investigation of much 


less studied finite temperature properties of anisotropic 
triangular quantum spin systems. Previous investiga¬ 
tions used variants of Monte Carlo methods for the quan¬ 
tum case IZH3, also for the classical case [laiii] or ana¬ 
lytical methods for the low temperature regime |12] and 
high temperature series expansion m- In the present 
work we use the finite temperature Lanczos method 
(FTLM) based on the exact diagonalization (ED) of fi¬ 
nite lattice tiles for this purpose. The method is based 
on evaluating the partition function by averaging over 
random starting vectors in the ED procedure. Using pe¬ 
riodic boundary conditions and discrete symmetries we 
can go up to largest cluster size of TV = 28 sites for ther¬ 
modynamic averages. We also consider smaller clusters 
of size N = 16, 20, 24 by FTLM. Furthermore we present 
the full analytical solution of the spectrum for the small¬ 
est iV = 8 cluster as a reference point. However, we 
note that unlike zero temperature ED results for finite 
clusters the FTL method cannot be used for finite size 
scaling analysis of thermodynamic quantities. Nonethe¬ 
less in the special case of the ID chain system (Ji = 0) 
we can compare to exact results which are in excellent 
agreement with the N = 28 cluster results. This indi¬ 
cates that our ETLM results are trustworthy to use as a 
representation for the thermodynamic limit. 

Our main emphasis is the calculation of peak position 
T/nax and height Cv{T,na.x) and xi^ms^x) in the temper¬ 
ature dependence of specific heat and susceptibility, re¬ 
spectively. We investigate its systematic variation with 
anisotropy control parameter J 2 /J 1 or cj) = tan“ ^ (J 2 / ), 

in particular when it is tuned between the simple special 
cases mentioned above. This information is of great im¬ 
portance for the analysis of thermodynamic data of trian¬ 
gular magnets. The Tniax(</>) dependence turns out to be 
considerably more asymmetric than for the related square 
lattice model [1]. We show that in certain ranges of the 
control parameter an indication of second maximum or 
shoulder in these quantities appears. Eurthermore as an 
example we discuss the analysis of the full temperature 
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FIG. 1. Schematic view of anisotropic triangular exchange 
models. For Ji = 0 (</> = ±'7r/2) this reduces to ID J 2 -spin 
chains (||), for J 2 = 0 (</) = 0) to the square lattice Ji- HAF 
(□) and for Ji = J 2 {<j) = rr/d) to the isotropic triangular 
system (A). 
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FIG. 2. Different exchange parametrizations used for the 
anisotropic triangular Heisenberg model and their dependence 
on the anisotropy parameter 0. The left and right branch for 
J 2 /J 1 corresponds to J 2 < 0 or J 2 > 0, respectively. 


dependence of x(T') for the most prominent triangular 
quantum spin systems CS 2 CUCI 4 and the isostructural 
Cs 2 CuBr 4 using FTLM results. We demonstrate that the 
derived exchange model is in excellent agreement with 
the one obtained from direct spectroscopic methods like 
inelastic neutron scattering (INS, CS 2 CUCI 4 only) spin 
wave results and electron spin resonance (ESR) experi¬ 
ments, both in high fields. 

In Sec. [n] we define the model and its parametriza- 
tion. The reference of the full solution for the eight-site 
cluster is discussed in Sec. mil The FTL method and its 
technical implementation are discussed to some extent in 
Sec. IV Our main results on specific heat and susceptibil¬ 
ity are presented in Secs. [V| and |VI[ respectively and the 
explicit comparison to CS 2 CUCI 4 and Cs 2 CuBr 4 is dis¬ 
cussed in Sec. |VII| Finally Sec. VIII gives the conclusion 
and outlook. 


II. ANISOTROPIC TRIANGULAR EXCHANGE 
MODEL AND ITS PARAMETRIZATION 

The anisotropic Ji — J 2 exchange model on the trian¬ 
gular lattice (Fig. is given by 

H = (I) 

{iJ) 

with 

T _ f If A 2 ^ 

\ J 2 if R, = R. ± e, ’ 

where ex,ey are unit vectors along cartesian directions. 
Figure may also be interpreted as tilted square lat¬ 
tice model with one diagonal J 2 bond cut out [ 2 ]. 
For the anisotropic triangular-lattice model, different 
parametrizations of the exchange energies are custom¬ 
ary. As introduced in Ref. [1] and used subsequently we 
prefer the polar representation defined by 

Ji = Jc cos (j), J 2 = Jc sin 4>, (3) 


Jc = + = tan 1 

where Jc gives the overall energy scale and </> is the 
anisotropy control parameter. This parametrization al¬ 
lows for an easy interpolation between important geo¬ 
metrical limiting cases, namely the square-lattice Neel 
antiferromagnet with J2 = 0 (((> = 0), the isotropic 120° 
triangular antiferromagnet with J 2 = Ji (^ = T^/d), the 
antiferromagnetic chain with Ji = 0 (^ = 7r/2), and their 
ferromagnetic counterparts. However, there are alterna¬ 
tive possibilities in the literature which we briefly men¬ 
tion here for ease of comparison (Fig. [^. Many authors 
regard the model as an extension to the one-dimensional 
spin chain, therefore they use the exchange J 2 along these 
chains as the overall energy unit and parametrize their re¬ 
sults in terms of a := J 1 /J 2 with Ji being the interchain 
coupling. From a square-lattice perspective in turn, the 
exchange parameter Ji appears to be the natural energy 
unit, and correspondingly 1/a := J 2 /J 1 is used as well. 

Both model parameters a and 1/a however are prob¬ 
lematic when trying to describe the full phase diagram, 
in particular the interpolation between the square-lattice 
antiferromagnet (a —>■ 00 ) and the one-dimensional chain 
(1/a —> 00 ). To overcome this, the function / := 
J 2 /( Ji + J 2 ) = 1/(1 + a) has been introduced in Ref. [H] , 
which remains finite in both limits and between these. 

However the latter parametrization also does not cover 
the full phase diagram unambiguously, which is why we 
have introduced the universal energy scale Jc and the 
anisotropy angle (j). In the square-lattice case, this has 
the additional advantage that Jc, apart from a global 
factor v^, denotes the magnetocaloric energy scale Jmc 
as well [1]. To facilitate the comparison of our results 
with the literature. Fig. displays the dependence of 
the quantities introduced here on the anisotropy control 
parameter (/. 

The nearest neighbor (n.n.) HAF model on the tri¬ 
angular lattice is the generic ‘geometrically frustrated’ 
spin system where the exchange energy of n.n. bonds 
cannot be minimized simultaneously for all bonds (see 
























3 


□ A I I 



0/7T 


FIG. 3. Degree of frustration in per cent in the triangular 
HAF as function of exchange anisotropy parameter. It van¬ 
ishes identical in the unfrustrated regime (p < 0 (J 2 < 0). The 
dotted lines are the classical phase boundaries [^. 


upper corner of Fig. [^. It is worthwhile to quantify this 
intuitive notion of ‘frustration’. A measure for it is the 
total loss of exchange energy due to frustration relative to 
the exchange energy without it. For the basic triangular 
(three-site) plaquette in Fig. the degree of frustration 
is then given by 


«(<(’) := 1 - 


Ea 

Et + Ed 


(4) 


Here Ea is the ground state energy of the frustrated 
triangle. Et := Ea{J 2 = 0) and Ed := Ea{Ji = 0) are 
the ground state energies of its unfrustrated trimer and 
dimer parts, respectively with 
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Ea{4') ■= min ( + -J 2 , -Ji + -J 2 


(5) 
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FIG. 4. Sketch of the tile 8:2-2. The edge vectors are ai and 
a 2 , the numbers label the individual sites. 


ing an averaging procedure over the low energy spec¬ 
trum only. To obtain a better intuition for the trian¬ 
gular Heisenberg model it is useful to calculate also the 
full spectrum by direct diagonalization of small clusters. 
We will demonstrate that it develops characteristic signa¬ 
tures at the classical phase boundaries and special sym¬ 
metry points as function of control parameter using the 
eight-site cluster. 


A. Spectrum and classical phase diagram 


As discussed in Ref. [3] for the square-lattice model, 
we can express the spectrum of the triangular eight-site 
cluster on tile 8:2-2 with periodic boundary conditions as 
a sum over Hamiltonians on complete graphs. We define 


Y. s,. 




i&C 


( 6 ) 


where the minimum is taken from the three different en¬ 
ergy eigenvalues of the triangle. Using these expressions 
in Eq. the degree of frustration (in per cent) in the tri¬ 
angular HAF as function of control parameter is shown 
in Fig. There the dotted lines indicate the classical 
boundaries between FM, AF and spiral phases as dis¬ 
cussed in Ref. [5]. «(</>) vanishes identically in the un¬ 
frustrated (J 2 < 0 or (/) < 0) case. In the frustrated 
regime (J 2 > 0 or ^ > 0) it achieves its maximum val¬ 
ues of K = 4/7 corresponding to 57% frustration at the 
isotropic point (A) and the Spiral/FM phase boundary 
while it vanishes for the ID chain case (||) where Ji = 0. 
We note that the reduction in the ordered ground state 
moment in the frustrated regime does not directly follow 
the degree of frustration but is the result of the subtle 
interplay of the latter with quantum fluctuations [5] • 


III. PRECURSOR: CLASSICAL PHASE 
DIAGRAM AND EIGHT-SITE FULL SOLUTION 

Within the finite temperature Lanczos method the 
thermodynamic quantities will be computed directly us- 


where C denotes an ordered list of lattice sites. In our 
case it has either two elements (a bond), four elements (a 
square or tetrahedron), or eight elements (a cube). We 
then can write 


Hs = Ji 


^CG 

«-{ 1 . 2 . 3 , 4 . 5 , 6 . 7 . 8 } 


(^CG 


5} 


+ H 


GG 

{ 2 , 6 . 7 , 8 } 


2J2 


n 


CG 

{ 1 , 5 } 


n 


CG 

{ 3 . 4 } 


-n 


CG 

{ 2 . 8 } 


n 


CG 

{ 6 . 7 } 


(7) 


for the Hamiltonian. The tile and its labelling is illus¬ 
trated in Fig. 1^ From Eq. Q with S = 1/2 being 
the maximum expectation value for each local spin S^, 
we get the corresponding list of eigenvalues by hierar¬ 
chically constructing all possible multiple-spin configura¬ 
tions starting with the basic two-spin singlets and triplets 
on the pairs {1,5}, {3,4}, {2,8}, and {6,7}: Erom all 
pairs of two-spin states, we can construct the four-spin 
states with total spin S = 0,1,2, each possible pair of 
these four-spin states in turn can be combined to an 
eight-spin state with total spin S' = 0,1, 2, 3,4 and spin 
degeneracy (2S -|- 1) . Because the total spin may be 
composed in several ways by the spins of sub-clusters 
there exist additional degeneracies. In this way, we can 
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TABLE I. All energy levels of the eight-site cluster. Sis is 
the total cluster spin and the next six columns give the sub¬ 
cluster spins. The last column denotes the additional degen¬ 
eracy of levels on top of the ( 25 'i _8 + l)-fold spin degeneracy. 
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construct a total of {j^/ 2 ) ~ states for the eight-site 
cluster. 

Table |T] displays the complete set of eigenvalues, to¬ 
gether with their (additional) degeneracies, total spins, 
and total spins on the sub-lattices. Of particular in¬ 
terest are those states corresponding to classical ground 
states (marked by an asterisk in the last column): the 
columnar antiferromagnet (first state in the table, en¬ 
ergy — 6 J 2 ), the Neel antiferromagnet (seventh state, en¬ 
ergy 2 J 2 — 6 J 1 ), and the ferromagnet (last state, energy 
2(2Ji -|- J 2 ). These three states replace each other as 
ground states as a function of the anisotropy ratio (/>, see 
Fig.i Fig. [^displays the full spectrum of the tile 8:2-2 
as a function of the anisotropy angle. The thicknesses 
of the lines indicate the degeneracies of the correspond¬ 
ing states as listed in Table Solid lines denote singlet 



FIG. 5 . Energy spectrum of the eight-site cluster, classihed 
according to total spin S. Solid lines: S' = 0 , short-dashed 
lines: S = 1 , dotted lines: S = 2 , dash-dotted lines: S = 3 , 
long-dashed lines: S = 4 . The line thickness indicates the 
degeneracy of the corresponding energy level, see Table [I] 
The thin vertical lines denote the classical phase boundaries 
between ferromagnet (EM), Neel antiferromagnet (AF), and 
spiral structure. The dotted vertical lines denote the three 
special cases square-lattice Neel antiferromagnet (Ji > 0 , 
J2 = 0 ), isotropic triangular lattice (J2 = Ji > 0 ), and one¬ 
dimensional antiferromagnetic chains (Ji = 0 , J2 > 0 ). 


states, the long-dashed line denotes the ferromagnet, for 
further line types see the figure caption. 

We overlay this spectrum to the classical phase dia¬ 
gram of the model discussed in detail in Ref. [2] . The clas¬ 
sical FM, AF and spiral phases are separated by the thin 
vertical lines. The thin dotted lines indicate as special 
cases the square-lattice Neel antiferromagnet (□, (j) = Q 
or Ji > 0 , J 2 = 0 ), isotropic antiferromagnetic triangular 
lattice (A, (j) = 7 r /4 or J 2 = Ji > 0) corresponding to 
the non-collinear 120 ° commensurate spiral structure and 
one-dimensional antiferromagnetic chains (||, </) = 7r/2 or 
Ji = 0, J 2 > 0). 

Characteristic for the spectrum at these six special 
points is the high number of degenerate excited states 
with different total spin. At the borders of the classical 
ferromagnetic phase also the ground state changes from 
the fully polarized state to a singlet. However the clas¬ 
sical phase boundary between the antiferromagnet and 
the spiral phase does not correspond to an equivalent 
change of (in this case nonmagnetic) ground states. This 
is due to the fact that with our eight-site tile we cannot 
approximate any incommensurate state at all, therefore 
the “best” approximation to the “true” ground state re¬ 
mains the Neel state for 1/2 < J 2 /J 1 < 3/4 (equivalent 
to 0.148 < (//tt < 0.205), and the columnar antiferro¬ 
magnetic state for larger values of J 2 /Ji. 

Similar to the square-lattice case |4], we observe level 
crossings near the classical phase boundary between AF 
and spiral phase. For the square lattice, we took this as 
an indication for the appearance of a nonmagnetic phase. 
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coinciding with the suppression of the ordered moment in 
a region around the classical transition observed within 
linear spin-wave theory [5]. Qualitatively the same hap¬ 
pens here, however this analogy should not be taken too 
far, because no indication of that kind exists for the clas¬ 
sically invisible crossover to one-dimensional chains at 
(j) = 7r/2, which is in reality surrounded by a large non¬ 
magnetic region. And empirically we know that linear 
spin-wave theory rather underestimates the size of these 
nonmagnetic regions. At the crossover from the spiral to 
the ferromagnetic phase for ^/tt « 0.852, a similar type 
of pattern of level crossings exists. Linear spin-wave the¬ 
ory gives inconsistent results in this case. 


B. Thermodynamics of the eight-site cluster 


Given the eigenvalue Ej with degeneracy dj and 
total spin Sj for each state I listed in Table 
we can easily evaluate the partition function Z = 
di (25'/ -I- 1) to calculate the magnetic suscep¬ 

tibility X s-ncl the specific heat cy according to 


X = 


Cy = 




N 


Y,di{2Si + l)E 


2^-l3Ei 


di{2Si + l)A/e 


-PEi 


(9) 


where /3 = l/lk^T) with kB being the Boltzmann con¬ 
stant. Here and in the following, we express the mo¬ 
lar susceptibility in units of NbHo where Nb 

is the Loschmid constant, /ig the magnetic permeabil¬ 
ity constant, g the gyromagnetic ratio, and /tb the Bohr 
magneton. For the dimension of the specific heat cy, we 
use the universal gas constant i? = Al/cb- 


tocaloric energy scale Jmc, which are dehned through 

= + ( 10 ) 

n 

= + (11) 

n 

The sums run over all bonds n connecting an arbitrary 
but fixed site i with its (not necessarily nearest) neigh¬ 
bors at sites {i + n}. We then have [IH] 

X = (1 - pe ), (12) 

cv = l[S{S + lfp^Jl,. (13) 

In principle we should be able to determine the exchange 
parameters Ji and J 2 already from high-temperature 
fits of the experimental results to the expressions above. 
However having fixed and 0 determines 2 J 1 + J 2 and 
I Ji — 4 J 2 I but leaves the sign of the latter undetermined. 
When expressed with the parameters Jc and p, this is 
equivalent to the fact that there are, apart from special 
cases, always two possible values (j)± for the anisotropy 
parameter. These two values (j)±, however, can lie in two 
different thermodynamic phases with completely differ¬ 
ent properties. This ambiguity is similar to the one in the 
square lattice Ji — J 2 model and its implications there 
were discussed in Refs. l4land[20l 

Although the coefficients of the high-temperature ex¬ 
pansions for x(T) and cv(T), being polynomial functions 
of Ji and J 2 , are known up to at least eighth order [TB] . 
this situation essentially will not change by including 
further higher-order terms in a high-temperature expan¬ 
sion [20) . and it remains difhcult to determine Ji and J 2 
unambiguously solely from fits to the high temperature 
dependence of x and cy. One powerful further diagnos¬ 
tic is the investigation of saturation fields um provided 
that they are in an accessible range. 


IV. GENERAL REMARKS ON FINITE 
TEMPERATURE METHODS 

The finite temperature properties of spin systems 
can be treated with analytical high temperature expan¬ 
sions [nuMm or with the numerical FTL method m 
which will be employed in this work. As a reference we 
first discuss briefly the single first order term of the ex¬ 
pansion method. 


A. High temperature approximation 

The high-temperature behavior of the susceptibility x 
and the specific heat cy to quadratic order in /3 is de¬ 
termined by the Curie-Weiss energy 0 and the magne- 


B. Finite-temperature Lanczos method 

To overcome this ambiguity, we use the finite- 
temperature Lanczos method [iiiiniiii] to evaluate the 
thermodynamic functions directly and compare specihc 
heat and susceptibility temperature dependence over the 
whole temperature range above the finite size gap re¬ 
gion. The method is based on the evaluation of thermo¬ 
dynamic traces using the eigenvalues and many-particle 
wave functions determined by numerical exact diagonal- 
ization of the Hamiltonian matrix on finite tiles. After 
mapping the Hamiltonian onto a sparse matrix, we use 
the iterative Lanczos algorithm [53] to generate the first 
few (between 1 and 100) extremal eigenvalues and the 
corresponding wave functions. Due to the Boltzmann 
weight, these are the most important eigenvalues con¬ 
tributing to the partition function. 
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FIG. 6. Tiles used in the finite-temperature Lanczos calculations. From left to right: 16:4-0, 20:2-4, 24:4-0, 28:2-4. The 
numbers label the lattice sites, ai and a 2 are the edge vectors. 


We classify our wave functions according to the ex¬ 
pectation value of the z component 
the total spin SI, the crystal momenta k and the point 
group symmetries of the tile. This brings the Hamilto¬ 
nian matrix into block diagonal form, allowing us to go 
to tile sizes up to = 28 in our finite-temperature cal¬ 
culations. In this way we can evaluate thermodynamic 
traces on industry standard computer hardware. 

The tiles we are using for the diagonalization proce¬ 
dures are illustrated in Fig. ^ We apply the same la¬ 
belling as described in Ref. [23] for the square-lattice case, 
adapted to the triangular lattice. There are two impor¬ 
tant differences: At finite temperatures we cannot easily 
perform any kind of finite-size scaling analysis of the par¬ 
tition function, therefore we directly take the results of 
the different lattice tilings. In addition there exists at 
least one phase corresponding to a classical phase with 
an incommensurate ordering vector. Due to the finite¬ 
ness of our k space grid, we cannot model this closely. 
However for thermodynamic properties like susceptibil¬ 
ity and heat capacity all thermally populated states con¬ 
tribute, and the exact modeling of the ground state as a 
function of our parameters is of secondary importance. 
For zero temperature ED approach the introduction of 
twisted boundary conditions may provide a way to cir¬ 
cumvent this problem |25j . 

We attempt to determine the thermal expectation 
value of an arbitrary static operator A by the fundamen¬ 
tal traces over the statistical operator, 

= (14) 

n—1 

■2^ = E ’ (1®) 

n—1 

where 2 is the partition function and Wt is the di¬ 
mension of the Hilbert space spanned by the basis 
{|n) : n = 1... Wt} (A^st ~ 2.7 • 10® for N = 28). In each 


symmetry-invariant subspace of the full Hilbert space, 
the Lanczos algorithm in principle is a sophisticated iter¬ 
ative basis change transforming the original Hamiltonian 
77 to an equivalent one for an open one-dimensional chain 
(not ring) problem with complex on-site potentials and 
nearest-neighbor interactions. Each iteration step corre¬ 
sponds to adding an additional site to this chain. After 
M steps, the resulting equivalent tridiagonal Hamilto¬ 
nian matrix T-Lm can be diagonalized easily to get the 
eigenvalues {tj : j = 1... M} and normalized wave func¬ 
tions {l^/’j) : j = 1... M} such that we can, in principle, 
evaluate Eq. (141. The moduli of the additional matrix 


elements of "Hm rapidly decrease with increasing number 
of iterations M, a property which we use as a conver¬ 
gence criterion. Eurthermore the difference in ground- 
state energy for iteration M and M — 1 is used as a 
second convergence criterion. Typical values such that 
machine precision is reached for the ground state energy 
are 10 < M < 100, which is vanishingly small compared 
to the original dimension Nst = O (lO®) of the Hilbert 
space. To sample an as large as possible part of the 
Hilbert space, we start A^r = 0(100) iterations with dif¬ 
ferent random wave functions or starting vectors |r), such 
that eventually we use no more than O (lO^) eigenvalues 
and wave functions per Hamiltonian block. It may be 
shown m that this procedure yields an asymptotically 
exact result. In summary, the thermal expectation value 
of A is approximated by 


, ly T3 ivjTD 

1 ATS R R 


R- r—1 j—0 




N, 


r=l j=0 


(17) 


where the index s denotes the summation over all symme¬ 
try sectors of the Hilbert space with dimension A))). If the 
operator A is a conserved quantity with [77, A] = 0, we 


can replace A by its quantum numbers AJ’®, and Eq. (16) 





7 


further simplifies to 


AT” 


r=l j=0 


The general definitions of the volume magnetic suscepti¬ 
bility and the heat capacity are 


d‘‘F 




(19) 


where H = Ba /Ho is the external magnetic field defining 
the z direction and F = (—l//?)lnZ is the canonical 
free energy. Written in terms of the expectation values 
defined above, we get the cumulants 




( 20 ) 

( 21 ) 


for the molar susceptibility and the specific heat of a tile 
with N sites, respectively. In our case, both and triv¬ 
ially H commute with H, allowing us to use Eq. (18) to 
determine the thermal expectation values. Furthermore 
spontaneous magnetic order is absent because our tiles 
are finite and we do not include an external magnetic 
field, therefore we can safely set {^z)p = 0. 

As mentioned before, in the comparison of calculated 
values and experimental results two strategies are possi¬ 
ble. One can identify the maximum position and value 
of thermodynamic quantities or perform a fit over the 
whole temperature range to extract the exchange model 
parameters. We will discuss both approaches. 


V. HEAT CAPACITY 

First we discuss the heat capacity which requires only 
the cluster eigenvalues for its evaluation. A compari¬ 
son of the theoretical temperature dependence to experi¬ 
ments is, however, not straightforward due to the lattice 
contribution to the heat capacity [3]. 


A. Temperature dependence 


To demonstrate characteristic features. Figs. [7] and 
show the temperature dependence of the specific heat 


cy(T) according to Eq. (21) for a selected range of 
anisotropy parameters (f> in the antiferromagnetic (Fig.[^ 
and the spiral phase (Fig. |^. The FTLM calculations 
are trustworthy only down to a temperature range of 
T « Jc/{NkB)- For lower temperatures they are dom¬ 
inated by the artificial finite size gaps. This excluded 
region is indicated by grey bars in the figures. 

For the AF phase with values —0.4 < 4>/Tr < —0.21 
a single peak indicated by dots that shifts continuously 


to higher values with increasing (j) is observed. This is 
qualitatively clear already from the eight-site spectrum 
(Fig§ which shows an increasing average excitation gap 
from the ground state in that range of (/). Furthermore 
a plateau (or very flat second maximum) at the lowest 
temperatures close to the finite size gap region is visible. 

For the spiral phase the values 0.25 < (/'/tt < 0.44, 
correspond to the region between the isotropic triangular 
lattice with J2/J1 = 1 and deep inside the disordered 
phase with quasi-one-dimensional chains {J2/J1 ~ 5). 
The tile 28:2-4 was used for the numerical calculations in 
both cases. 

In Fig. ^ the characteristic maxima of cv(T), indi¬ 
cated by the small black dots in the figure, have positions 
which fall into two different temperature ranges: For the 
isotropic triangular case with (jj/ir = 0.25, the maximum 
sits at Ta ~ 0.18Jc/fcB, reducing in magnitude and shift¬ 
ing slightly to lower temperatures Ta « 0.15Jc/fcB upon 
increasing </> from its isotropic value to (jj/ir « 0.31 or 
J2/J1 ~ 1.5. At this point, a second maximum devel¬ 
ops starting at Ty « 0.32Jc/fcB, shifting to higher tem¬ 
peratures as 4> is increased towards the disordered re¬ 
gion. Furthermore, the maximum Ta characteristic for 
the triangular lattice rapidly disappears with increasing 
4 >. Only the Cv{T) curves for ^/tt = 0.31 and ^/tt = 0.32 
show both maxima simultaneously. At high temperatures 
T 3> Tc/Zcb, all the curves show the (1/T)^ temperature 
dependence expected from Eq. (13). 


B. Exchange anisotropy dependence of peak 
position and value 


Fig. [I shows a compilation of the maximum positions 
and values of cv[T) for the five different tiles we have 
used as a function of the anisotropy parameter (j). The 
solid lines denote the exact solution for the tile 8:2-2, 
see Eq. and Table [I] The black dots denote the max¬ 
ima for the tile 16:4-0, the open triangles for the tile 
20:2-4, the solid diamonds for the tile 24:4-0, the open 
squares for the tile 28:2-4. Furthermore, the exact re¬ 
sults taken from Ref. [26] for the one-dimensional chain, 
cv(Tinax)/(A^L*B) = 0.35 and /cBT^ax/^c = 0.481, are 
marked with the white circles and serve as a gauge to 
judge how close a cluster of given size approaches the 
thermodynamic limit. 

The eight-site cluster, introduced as an illustration 
of the model and its overall spectrum, clearly is not 
very useful quantitatively to discuss the heat capacity 
in the thermodynamic limit. Only the overall behavior 
of the maximum temperature is qualitatively similar to 
our findings for the larger tilings, however the two min¬ 
ima in the spiral phase are shifted towards the classical 
phase borders. The values cv{Tjnax) for the eight-site 
cluster are monotonically increasing in the antiferromag¬ 
netic phase, followed by a minimum at (p/n « 0.1 and 
a maximum at the isotropic point, (pjTr = 0.25. On the 
ferromagnetic-Ji side for 1/2 < p/n < 1, cv(Tmax) in- 
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FIG. 7. Temperature dependence of the specific heat cv{T) of the anisotropic triangular lattice according to Eq. ( |21[ | for 
anisotropy parameters 4> in the antiferromagnetic phase ranging between tp/ix = —0.4 and ^/tt = —0.21, see legend. We used 
tile 28:2-4 for the numerical evaluation of Eq. ( |21[ ), the grey-shaded area at low temperatures illustrates the finite-size gap of 
order 0{Jc/N). The characteristic maxima of cv(T) are indicated by the small black dots. 
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FIG. 8. Temperature dependence of the specific heat cv{T) of the anisotropic triangular lattice according to Eq. (21 1 for 
anisotropy parameters (f) in the spiral phase ranging between (j>/TT = 0.25 (isotropic triangular lattice) and c^/tt = 0.44 (crossover 
to one-dimensional chains), see legend. We used tile 28:2-4 for the numerical evaluation of Eq. ( |21[ ), the grey-shaded area at 
low temperatures illustrates the finite-size gap of order 0{Jc/N). The characteristic maxima of Cv{T) are indicated by the 
small black dots. 


creases again to a second maximum, followed by a min¬ 
imum at the crossover to the ferromagnet at J 2 /\Ji\ = 
1/2. At the one-dimensional point Ji = 0, naturally both 
Tmax and cviT^a^) differ strongly from the exact values. 

Also the 16-site tiling shows similar features. For the 
larger clusters of size N = 20 and more, both max¬ 
imum temperature positions and maximum values of 
the specific heat clearly show a double-peak structure 
as function of </. The positions Tmax decrease with in¬ 
creasing cluster size, apart from the regions around the 
isotropic point and in the spiral phase near the ferro¬ 
magnet. At the one-dimensional point, agreement with 
the infinite-chain result [26] (white circles) is achieved 
already with TV = 20 tiling. And therefore results for 
the N = 28 tile represent well the thermodynamic limit 
as far as peak position and height is concerned. De¬ 


noting the largest value in each sector of Fig. m by 
'T'max we observe that T,/,,,,,(AF)/r/,,,^(spiral) = 1.58. 
This asymmetry in Tiiiax(</') is considerably larger than 
the corresponding one in the square lattice J 1 -J 2 model 
where T,^ax(^F)/^max(CAF) = 1.26. (The columnar AF 
(CAF) phase replaces the spiral phase in this model.) 
Similarly for the asymmetry ratios of peak values in the 
triangular case we get cv[T^^J{AF)/cv[T^^J{spiTal) = 
1.51 much larger than cv[T^i^.^]{AF)/cv[T^^^]{CAF) = 
1.10 in the square lattice. Inside the AF and spiral sectors 
the lowest Tmax is reached close to the isotropic triangu¬ 
lar point with T^ax ~ O.lSJc/fce because there exchange 
frustration is most pronounced (Fig. [^. This leads to a 
high density of low energy excitations (c.f. Fig. and 
therefore a low 

Furthermore, in parts of the antiferromagnetic as well 
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FIG. 9. Maximum positions (top) and values (bottom) of the 
specific heat cv (T) as a function of the anisotropy parameter 
<j>. Solid lines: exact solution for tile 8:2-2, squares: tile 16:4- 
0, triangles: tile 20:2-4, diamonds: tile 24:4-0, filled circles: 
tile 28:2-4. The white rings mark the corresponding exact 
values for the one-dimensional chain |26| . 


FIG. 10. Maximum positions (top) and values (bottom) of 
the static susceptibility x{T) as a function of the anisotropy 
parameter (f). Solid lines: exact solution for tile 8:2-2, squares: 
tile 16:4-0, triangles: tile 20:2-4, diamonds: tile 24:4-0, filled 
circles: tile 28:2-4. The white rings mark the corresponding 
exact values for the one-dimensional chain | 26| . 


as the spiral phase, a second maximum cy(rmax 2 ) at 
very low temperatures appears. T'niax 2 depends roughly 
linearly on (j) while cy (Tniax 2 ) remains constant and 
small, see also Fig. With the 28-site tiling being the 
largest possible, we cannot judge whether this second-low 
temperature maximum in the antiferromagnetic phase 
merely is a finite-size effect. This is different from the sit¬ 
uation in the spiral phase, where the partially observed 
second maximum is at temperatures rmax 2 far larger than 
the finite-size gap, see also Fig. 

Similar to the behavior around the isotropic point, at 
the crossover to the ferromagnetic region a second max¬ 
imum gradually appears in cy(r) which then evolves to 
the “ferromagnetic” maximum while the “spiral” max¬ 
imum disappears. This also happens at comparatively 
low temperatures with an irregular behavior of the max¬ 
imum values as function of 4>. We attribute this irreg¬ 
ularity to the fact that in particular near the borders 
of the spiral phase the ground state and possible low- 


lying excited states have incommensurate ordering vec¬ 
tors Q = Qaf + i 5Q and Qfm + i 5Q respectively with 
|5Q| 1. These are not contained in the coarse grid of 

crystal momenta of our finite tiles. 


VI. SUSCEPTIBILITY 


Evaluation of the susceptibility with FTLM is more 
involved because it requires the calculation of magnetic 
moment matrix elements. On the other hand this quan¬ 
tity can be more easily compared to experimental results. 
Most known exchange parameters for spin systems are 
due to the application of this method. 
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FIG. 11. Temperature dependence of the static susceptibility x{T) of the anisotropic triangular lattice according to Eq. ( |20[ ) 
for anisotropy parameters (f) in t he s piral phase ranging between = 0.57 and (p/iv = 0.86, see legend. We used tile 28:2-4 for 
the numerical evaluation of Eq. (201, the gray-shaded area at low temperatures illustrates the hnite-size gap of order 0{Jc/N). 
The characteristic maxima of x{T) are indicated by the small black dots. 


A. Temperature dependence 


Similar to the specific heat, we have calculated the 
temperature dependence of the static magnetic suscepti¬ 


bility x(T) according to Eq. (20) with the tilings shown 
in Fig. Fig. [TT] shows the results in a selected range of 
the anisotropy parameter, 0.57 < (ji/ir < 0.86, calculated 
with tiling 28:2-4. This parameter range corresponds 
to the interpolation between ferromagnetically coupled 
(Ji < 0) quasi-one-dimensional antiferromagnetic chains 
with JijJ\ ~ —4.5 and the crossover to the ferromagnet 
at the classical boundary J2/J1 = —1/2. For the former, 
a very broad maximum is characteristic which gradually 
evolves into the T = 0 divergence of x(r) in the ferro¬ 
magnetic phase, which is the expected behavior. 


B. Exchange anisotropy dependence 

In the same way as for the specific heat, we follow 
the positions and values of the characteristic maxima of 
x(T) as a function of the anisotropy parameter, using the 
tiles displayed in Fig. again. This is shown in Fig. [T0| 
The solid lines denote the eight-site results according to 
Eq. For the same reasons as discussed for the heat 
capacity, strong deviations from the larger tilings occur 
in particular in the spiral phase. In contrast to cv(T), 
the maximum positions and values of x(^) already 
converged for tilings of size N = 20 and larger, apart from 
the crossover to the ferromagnetic phase, where a tile- 
dependence of the maximum temperatures clearly can 
be observed. 

Similar as for the specific heat the Ti„ax((/') depen¬ 
dence has a pronounced asymmetry in AF and spi¬ 
ral sectors. Denoting again the largest value in each 
sector by obtain for the triangular lattice 


’^max(f^F)/^max(spi'^al) = 1-58 which is much larger than 
the asymmetry value T,^ax(AF)/’^max(CAF) = l.Il for 
the square lattice case. Similar to specific heat behav¬ 
ior the T'i„ax(</) minimum inside AF and spiral sectors is 
reached around the most frustrated isotropic triangular 
position « O.SbJc/fcs- 

At the one-dimensional point, as for the heat capac¬ 
ity, the results from Ref. [5^, feeTmax/^/c = 0.641 and 
x(T’max)/(AfLMo(ff7‘B)^/-/c) = 0.147 are accurately repro¬ 
duced with our method. In general we think that in 
particular our results on the magnetic susceptibility can 
be used to accurately determine both exchange constants 
Ji and J 2 individually, which we show in the following 
section. 


VII. APPLICATION TO CsaCuCU AND 
Cs2CuBr4 


As a demonstration of the usefulness of the method, 
we apply our findings to the compounds CS2CUCI4 and 
Cs2CuBr4. Fig. [T^ shows the temperature dependence of 
the magnetic susceptibility. The dots and the open cir¬ 
cles denote the experimental data taken from Ref. 
the solid and dashed lines denote the corresponding fits 
with our FTLM data. In order to avoid misunderstand¬ 
ings, we have plotted the data and the curves in two 
different ways: The top plot in Fig. 12 displays x{T) in 


electromagnetic (CGS) units from where we have deter¬ 
mined the values for our model parameters reproduced 
in Table [TT] Using these fitted values for Jc, </, and g, 
we have plotted the same data again in dimensionless 
units in the bottom plot of Fig. The main effect 
of replacing Cl with Br is an increase of the overall en¬ 
ergy scale Jc by a factor 3.4, whereas the anisotropy an¬ 
gle (j) changes only by about 5 %. Thus only this large 
change in Jc is responsible for the decrease of the maxi- 
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TABLE II. Comparison of exchange parameters for CS2CUCI4 and Cs2CnBr4 as determined from thermodynamic FTLM-fit 
and direct spectroscopic INS (CS2CUCI4 only) and ESR methods aX H > Hsa.t- Exchange constant J2 corresponds to the 
crystallographic b direction and Ji to the zigzag bonds in the be plane. 


compound 

method 

Ji / meV 

J2 /meV 

Jc/meV 

J2/J1 


9 

Ref. 

CS2CUCI4 

FTLM 

0.11 

0.38 

0.40 

3.45 

0.41 

2.06 

this work 


INS 

0.128 

0.374 

0.40 

2.92 

0.40 

2.19 

m 


ESR 

0.12 

0.41 

0.43 

3.42 

0.41 

2.08 

m 

Cs2CuBr4 

FTLM 

0.5 

1.26 

1.35 

2.52 

0.38 

2.04 

this work 


ESR 

0.53 

1.28 

1.38 

2.44 

0.38 

2.09 

m 



kgTIJc 

FIG. 12. Temperature dependence of the magnetic suscepti¬ 
bility x{T) of CS2CUCI4 (dots and solid line) and Cs2CuBr4 
(open circles and dashed line). The top plot displays the ex¬ 
perimental data (symbols, taken from Ref. [^) together with 
the fits of our FTLM data (lines, fitted values see Table [II|). 
The bottom plot displays exactly the same data, this time in 
dimensionless units using Jc, (j), and g for the two compounds 
from Table ini 


mum in the broadening of the maximum, and the 

shift of the maximum towards higher temperatures. The 
anisotropy ratio and therefore the position of Cs2CuBr4 
in the phase diagram remains essentially the same as for 
CS2CUCI4, it is only slightly moved away from the quasi- 
one-dimensional regime. Historically the first measure¬ 
ment of x(T) of CS2CUCI4 yielded slightly different val¬ 
ues [30] , however the data analysis was performed over a 


limited temperature range assuming weakly coupled one¬ 
dimensional chains. 


As pointed out in Sec. V A the steep decrease of x{T) 
for T —)■ 0 in both cases is an artifact of the finiteness 
of the tiling used for the FTLM calculations. Assuming 
a finite-size gap A « Jc/N, the corresponding tempera¬ 
tures where we expect finite-size effects to dominate are 
Ta « 0.2K for CS2CUCI4 and Ta « 0.6K for Cs2CuBr4. 
Experimentally, the lowest temperature was Tmin « 2 K, 
well above the finite-size gaps. 

Our result from the FTLM fit to thermodynamic data 
are compared in Table |lT] to results from direct spectro¬ 
scopic methods: Inelastic neutron scattering (INS) [33] 
(CS2CUCI4 only) and electron spin resonance (ESR) [29] . 
both in fields above the saturation field [3] 


A^O-^sat — 


2S (Ji -b 2 J2)^ 


5MB 


2J2 


( 22 ) 


which is /J.oi?sat ~ 8.4 T for CS2CUCI4 and /io^^sat ~ 
SOT for Cs2CuBr4. Taking the spectroscopic data at 
high fields has the advantage that the ground state is 
fully polarized, corresponding to a single “all spins up” 
spinor product state |FM). The excitations on top of 
this are the N orthonormal single-particle excitations 
IV-i) = (l/v^)5'-|FM), i = 1...N with N = 0{Nl). 
Because [H, Hz] = 0, the Hamiltonian can not generate 
more than these one-spin-flip states, and its spectrum 
can be determined exactly by Fourier transform. 

The agreement between the different methods is almost 
perfect, for exchange parameters as well as g-factors. As 
already mentioned in Ref. [3| if CS2CUCI4 is interpreted 
as a purely 2D system this set of exchange parameters 
puts the compound very close to the quasi-ID spin liquid 
regime centered at </> = O.Stt. In reality, however, the 
magnetic order is stabilized below Tn = 0.62 K by a 
finite inter-plane coupling of the order J± « 0.017 meV 
along the crystallographic a direction |3T] . 

For Cs2CuBr4, we see deviations of the experimental 
data from our result at the lowest temperatures. These 
are not due to impurities [33133], but can be regarded as 
an indication for a tendency towards magnetic order in 
this compound. The overall energy scale Jc is more than 
three times larger than for CS2CUCI4, however the esti¬ 
mate for the anisotropy angle (j) is essentially the same. 
We therefore expect that Cs2CuBr4 orders magnetically 
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as well, possibly at a higher temperature than its Cl coun¬ 
terpart. 

The agreement of our thermodynamic analysis with 
direct spectroscopic results gives us confidence that the 
FTL method may be applied to the analysis of the whole 
Cs2CuCl4_a;Brx substitutional series [2ai32]. 


VIII. CONCLUSION AND OUTLOOK 


In this work we have shown that finite temperature 
Lanczos method for finite clusters is a versatile tool to 
investigate the little known finite temperature properties 
of frustrated triangular quantum magnets, in contrast to 
QMC approach which is not suitable in this case. The 
FTL method complements the analytical spin wave ap¬ 
proach which is useful only for the very low temperature 
regime and is a more straightforward alternative to the 
high temperate series expansion method. 

For this quantum spin model one can use a single con¬ 
trol parameter and tune the system through the phase di¬ 
agram where several special cases with very different frus¬ 
tration degree like frustrated Neel, unfrustrated HAF, 
frustrated spiral and isotropic triangular 120° phase as 
well as unfrustrated spin chains and frustrated FM may 
be realized. 

For the largest investigated cluster with N = 28 sites 
we obtain a trustworthy representation of the thermo¬ 
dynamic limit behavior since the exact known result for 
the spin chain case is produced very well and shows little 
difference to the iV = 24 cluster. As a main result we 
gave the systematic variation of peak position and peak 
height of specific heat and susceptibility as function of 
J2/J1. Both values are strongly suppressed for the most 
frustrated isotropic triangular magnet (Figs. and 10). 


Furthermore a surprisingly large asymmetry in partic¬ 
ular for the maximum position exists between the AF 
and spiral phase region. It is considerably larger than 
between the AF and CAF regions in the square lattice 
J1-J2 model. In addition we found that the simple single 
peak shape of cv{T) and x(T) may be distorted due to 
smaller side maxima or shoulders at lower temperature 
for some ranges of the control parameter. 

The most reliable method to extract the exchange pa¬ 
rameters is the fitting of x(T) over the whole tempera¬ 
ture range (above the finite size gap) using FTLM results. 
We have demonstrated this for two of the most typical 
2D triangular magnets, CS2CUCI4 and Cs2CuBr4. We 
have obtained excellent agreement with results from the 
INS investigation and with ESR results in the fully po¬ 
larized state. Our method has the additional advantage 
that it can easily be extended to the whole substitutional 
series Cs2CuCl4_a;Brx (0 < x < 4) [27] to extract the 
systematic variation of frustration anisotropy control pa¬ 
rameter (j) and energy scale Jc as a function of Br con¬ 
centration [33]. We note that our method can also be 
applied to magneto caloric measurements of the organic 
charge-transfer salts where localized S' = 1/2 magnetic 
moments are residing on a possibly distorted triangular 
lattice as well. This includes the Eta;Me4_a;Z[Pd(dmit)2]2 
(x = 0,1, 2) family of compounds [M] [3S| as well as k- 
(ET)2B(CN)4 and k-(ET) 2 Cu 2 (CN )3 |3Sj and k-(BEDT- 
TTF)2Cu[N(CN) 2]C1 [37j. 
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